Density matrix renormalization group algorithms for Y-junctions 
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Systems of Y-junctions are interesting both from a fundamental viewpoint and because of their 
potential use in nanoscale devices. These systems can be studied numerically with the density 
matrix renormalization group(DMRG), but existing algorithms are inefficient. Here, we introduce 
a much more efficient DMRG algorithm for Y-junction systems. As an example of the use of this 
method, we study S = i bound states in Heisenberg 5 = 1 junctions with two geometries, one 
where the junction consists of a single site, and the other where it consists of a triangle of three 
sites. 



Junctions are essential ingredients of existing and fu- 
ture electronic and spintronic devices. As progress to- 
wards smaller and smaller devices continues, eventually 
we will reach the atomic scale. Perhaps the smallest con- 
ceivable junctions are composed of intersecting spin or 
electron chains, with each site representing a single atom. 
Isolated chains, modeled by Heisenberg, Hubbard, and 
related Hamiltonians, have been extensively studied. Yet 
surprisingly little work has been done on the correspond- 
ing junctions of three or more half-chains, or legs. There 
has been recent work on junctions of quantum wires. For 
example, Oshikawa, et a/Q,Q studied the transport prop- 
erties of a junction composed of three quantum wires en- 
closing a magnetic flux, modeling the chains as spinlcss 
Tomonaga-Luttinger liquids, and finding a rich phase di- 
agram. Quantum Monte Carlo on junctions of fermion 
chains are subject to the minus sign problem, so little 
numerical work has been done. We are not aware of any 
studies which have focused on junctions of spin chains. 

Our present work has two main aspects. First, we 
develop a new density matrix renormalization group 
(DMRG) algorithm for studying junctions. Ordinary 
DMRG[3] is much less efficient for junctions compared to 
chains with open boundaries; our new method requires 
only slightly more computational effort. Second, as an 
application of this method we consider a key feature of 
a few of the simplest junctions, namely the presence of 
S = \ "spinon" bound states in two types of S = 1 
Heisenberg Y junctions. These bound states are the ana- 
logues of the spin-i end states found in open S = 1 



chainsjj, |a|- They exist at any S = 1 junction with an 
odd number of legs, regardless of the local interactions at 
the junction. Magnons coming to the junction can inter- 
act with the bound state via spin-flip scattering, leading 
to potentially interesting dynamical phenomena. We call 
these spinon bound states because a spin- i degree of free- 
dom is bound to the junction, but it is important to note 
that they are features of the ground state. A separate is- 
sue, which we do not address here, is whether the lowest 
excited state is below the Haldane gap, representing, say, 
a magnon bound to the junction. For example, a single 
5=2 impurity in a S = 1 chain has such a bound state 
below the gap only for sufficiently weak coupling to the 
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FIG. 1: Schematic view for the geometric configuration and 
algorithm for Y-junctions. The geometry of the Y-junctionis 
shown in (a), where the circles represent sites and the solid 
lines represent connections. The alternative geometry Ya 
junction is shown in (b). The original DMRG method is 
shown in (c), where one cuts the link between two sites and 
relinks it to the center of the remaining chain, making it topo- 
logically equivalent to a Y-junction. Verstraete and Cirac's 
tensor representation is shown in (d), where each open cir- 
cle represents a tensor and each solid dot represents a tensor 
index. An Affleck-Kenedey-Lieb-Tasaki valence bond picture 
of the Y-junction is shown in (e), where each solid dot rep- 



resents a S 



the line segments represent valence bonds, 



and where each open circle represents a projection operator 
to the S — 1 subspace for a site. 



impurity [g. 

We consider two types of three-leg junctions, which we 
will label Y and Ya, where we think of the A as a tri- 
angle. These are shown in Fig. QJa) and (b). DMRG 
is extremely efficient for one dimensional systems with 
short range interactions. It is standard to treat two- 
dimensional strips by mapping the strip into a chain with 
interactions over a distance of several lattice spacings. 
Similarly, a Y junction can be treated as a chain with 
one long range interaction, shown in Fig. W[c). The 
long-range interaction is equivalent to periodic bound- 



ary conditions in terms of the effort necessary in DMRG: 
a typical block of sites connects to the rest of the sys- 
tem by two links in both cases. If to states per block are 
required for a given accuracy for a single chain, one ex- 
pects roughly to 2 states to be needed for the Y junction 
(or a periodic chain). This changes the computational 
effort from 0(Nm 3 ) to 0(Nm 6 ), where N is the number 
of sites. 

As discussed below, our new algorithm scales as 
0(Nm 3 + m ). Typically to and iV are similar in size 
(say ~ 10 2 ), so the cost is close to that of a single chain, 
and a large improvement over the method with a long- 
range interaction. Our algorithm is somewhat related to 
Otsuka's early DMRG treatment of the XXZ model on 
a Bethe lattice |j|. It is closely related to Shi, Duan, and 
Vidal's very recent treatment |8| of general tree-tensor 
networks using Vidal's time-evolving block decimation 
method, which is closely related to DMRG. One can also 
view it as perhaps the simplest step towards an imple- 
mentation of Verstraete and Cirac's much more compli- 
cated tensor approach [SUlfll, which holds great promise 
for true two dimensional systems. 

Consider the matrix product representation for the 
wavefunction in standard DMRG|ll| 



li> 



lk> 



|f) = A[s 1 }A[s 2 }A[s3}-A[s n }\s 1 S2S 3 ...s r , 



(1) 



where Sj labels the states of site i, and A[si] is a matrix 
except for the first and last sites, where it is a vector. One 
can view each matrix as a two-dimensional tensor, with 
indices connecting to the left and right. In Verstraete 
and Cirac's approach for a square lattice, one generalizes 
to four-dimensional tensors, with indices connecting to 
each nearest neighbor. In our approach, only one ten- 
sor is needed, at the junction site, and it has only three 
indices. All sites on the legs are ordinary matrices, as 
shown in Fig. ^d). Unlike in the more general tensor 
approach, the states can be easily described in terms of 
an orthogonal basis. 

The key to our approach is a new type of step for the 
junction. In this "junction step" the system is divided 
into three blocks instead of two (one for each leg), and 
the calculation time is higher, <3(to 4 ), where to is the 
number of states kept per block, instead of 0(m 3 ). The 
rest of a sweep consists of steps involving two blocks, 
where one of the blocks contains two legs plus part of the 
third. The sweep moves out to the end of one leg, and 
then back to the center, after which the junction step is 
repeated, and then the sweep goes out another leg, etc., 
until all legs have been treated. 

We now describe this method in more detail. The sys- 
tem is built up in a warmup sweep. An ordinary DMRG 
chain initialization sweep is first carried out on one leg of 
the Y-junction. For the environment block one can use 
the next several sites of the leg, or artifically reflect the 
leg-the algorithm is not very sensitive to the warmup, 
and we keep only a small number of states to, storing the 
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FIG. 2: Configurations illustrating the new DMRG algo- 
rithm. The procedure starts from the center (a), where two 
legs are combined to form a new system block. The sweep 
next progress down a single leg, as in (b), reaches the end of 
the leg, (c), and comes back to the center. 



related operator matrices and transformation matrices. 
One continues until the entire leg is incorporated into 
the block. At the last few steps the environment may or 
may not be in the form of a junction — it does not make 
much difference. This block is then duplicated twice to 
form the other two legs. 

At this point, the finite-system sweeps begin. The first 
step is the junction step. For the Y-junction, where the 
junction consists of a single site, that site is incorporated 
into one leg. The wavefunction can be written as ipijk, 
where i, j, and k refer to the three legs/blocks. Ap- 
plying the operators appearing in the Hamiltonian H to 
ip requires a calculation time of O(to) 4 , in contrast to 
the O(m) 3 time on a chain. For example, to apply the 
term 5^,5?, connecting the first two legs, one first mul- 
tiplies the matrix 5? , into VVj'fc, storing the result, and 
then multiplies that result by Sf^ . Several such multi- 
plications by H are performed to make a few Lanczos or 
Davidson steps to find an approximate ground state ip. 

After finding ip we want to combine two legs (say j and 
k) into a single effective block, as shown in Fig. [3(a). In 
ordinary DMRG, we find the reduced density matrix for 
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FIG. 3: Error in the energy as a function of the number 
of states kept m for the new and the standard DMRG algo- 
rithms. The reference energy E here is the result of standard 
algorithm by keeping 3000 states. 



FIG. 4: Results for (Sz(l)) for a Y-junction as a function of 
site I along a chain, compared with an open ID chain. For 
the Y-junction, the first site shown is the center junction site. 
Here, the total spin of the system is S z — i. 



the system block, and diagonalize it, which in this case 
would result in a calculation time of 0(m e ). Instead, 
we perform a singular value decomposition (SVD) on ip, 
treating j and k together as a single combined index, 
and i as the other index of the matrix. The SVD has 
a calculation time of only 0(m 4 ), but is equivalent to 
the density matrix diagonalization. (Note, however, that 
the SVD is not as easily generalized to having more than 
one state targeted.) One obtains transformation matrices 
O a ,jk, where a are the new states of the combine blocks. 
The operators for the combined block A' are obtained as 
OAO T ; this can also be performed in 0(m 4 ) time. For 
example, to transform Sf^S^,, one first forms X a jr k = 

O a jkSjj,, then Y a ijik = S^ k ,O a >j>k', and then contracts 
XY over indices j'k. 

Next, as shown in Fig. HJb), we sweep down leg i. 
This sweep is an ordinary DMRG sweep, with calculation 
time 0(m 3 ) per site; the fact that the system block now 
includes two legs does not affect calculations in each step. 
We continue moving to the edge of the leg, Fig. Hfc), 
and then turn around and come back to the center. The 
junction step is performed and then the sweep moves out 
one of the other legs. The three legs need not be identical; 
the sweep treats each independently. In each full sweep, 
the center step is performed three times. Alternatively, 
if the Hamiltonian and the desired state are symmetric 
with respect to the interchange of legs, just before the 
junction step one can replace each of the legs with the 
most recently updated leg/block. 

For an 3 x N system, the total cpu time is 0(m A + 
Nm 3 ). In comparison, we consider the standard DMRG 
calculation with one long range interaction. In Fig. [3]wc 
show the error in the energy as a function of the number 



of states kept for the new versus the standard algorithm; 
the new algorithm requires vastly fewer states. If we 
assume the long range interaction mandates 0{m 2 ) states 
be kept for sufficient accuracy (for large N), the standard 
DMRG calculation would scale as 0(Nm 6 ). 

We now use the new algorithm to study S = | bound 
states in S = 1 Heisenberg junctions. A simple consid- 
eration of quantum numbers makes the presence of the 
bound states clear. An S = 1 open chain has S — \ states 
on the ends, and a finite correlation length £ = 6.03 |j|. 
The ends of the legs in a large junction system will also 
have S = \ states; the finite correlation length does 
not allow the junction to influence the ends. With an 
odd number of legs, a half-integer-spin state must form 
at the junction to make the total spin an integer. For 
a junction formed from very weakly coupled legs, the 
bound state comes from coupling three effective S = 4 
states at the junction. For more strongly coupled legs one 
can think in terms of the Afflcck-Kcnncdy-Lieb-Tasaki 
(AKLT) picture 12] where each S = 1 is composed of 
two S = i's, with one valence bond singlet connecting 
two S = |'s on each near-neighbor link; see Fig. 1(e). 
At the junction, there must be one S = | left over after 
the singlets are formed. 

In order to study the junction bound state, it is conve- 
nient to place a real S — \ spin on the end of each chain 
(far from the junction), which forms a singlet with the 
effective S = |, eliminating it as a low energy degree of 
freedom. For the Y-junction, the ground state is then a 
nondegenerate doublet corresponding to total S z = ±tj- 
In Fig. 0]we show (S z (l)) as a function of site / for a 
3 x 60 Y-junction. The result for the other two legs is 
identical. We see a strong similarity with the results for 
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FIG. 5: Results for (Sz(l)) for a 3 x 60 YA-junction, where 
site is part of the junction triangle. 



the end state of a single chain, except for a smaller mag- 
nitude of the oscillations. The decay of the oscillations is 
consistent with an exponential decay with the correlation 
length of a single chain|j|. Looking near the junction we 
find the initially surprising result that the value of {5 , JZ } | 
at the junction site (0.30) is less than that on the first 
site of a leg (0.34). If one thinks in terms of a mean- 
field/classical picture, with the spin fluctuations acting 
like thermal fluctuations, one would expect the opposite: 
the junction site is surrounded by three polarized spins, 
and would experience a larger polarizing field than the 
adjacent sites which have only two neighbors. However, 
if one utilizes the AKLT valence bond picture instead, 
one would predict a smaller magnitude at the junction 
site: the extra unpaired 5 = \ is energetically favored to 
be on an adjacent site, with one missing valence bond, 
as shown in Fig. ^e). An unpaired spin on the junction 
site occurs from fluctuations of the unpaired 5 = | from 
leg to leg. 

Now consider the other geometric configuration, the 
YA-junction. Fig. shows (S z (l)) for a YA-junctionwith 
size 3 x 60. One sees the presence of a bound state at the 
junction, as for the Y-junction. However, the symmetry 
under interchange of legs is broken by an exact ground 
state degeneracy. 

To understand this degeneracy, consider a three-site 
5=5 Heisenberg triangle. The S = \ triangle is clearly 
relevant in the limit that the exchange couplings connect- 
ing the three junction sites to each other is very small. 
In that case the S = -k 's are the end states of each leg. 
In the opposite limit, where the exchange between the 
5 = l's in the junction is large, the ground state of the 
5 = 1 triangle is nondegenerate (with a gap of precisely 



J), but then one expects a 5 = \ end state on each 
leg starting at the site adjacent to the junction, and one 
would expect these three S — ^'s to be weakly coupled 
antifcrromagnetically through the central triangle. 

If we square the total spin operator for a generic 
Heisenberg triangle, we obtain 



E 
1 



1 



i, and 5tot can be either \ or |. 
j J and consist of two 



(5 tot (5 tot +l)-5 1 (5 1 +l)-5 2 (5 2 +l)-5 3 (53+l)). 

(2) 
Here Si = 5 2 = 5 3 
The ground states have energy 

degenerate 5 tot = \ doublets. If one regards the triangle 
as a periodic chain, then the degeneracy corresponds to 
total momenta around the triangle of k = ±2ir/3. The 
5tot = | multiplet has k = and energy | J. The YA- 
junction has the same degeneracy as the ground states 
of the triangle; the numerical calculation converges to a 
arbitrary linear combination of these states. We have 
confirmed the degeneracy behaves as expected by tar- 
getting several states on a relatively small YA-junction 
system. 

In conclusion, we have presented a relatively simple 
DMRG algorithm for junction systems, with greatly im- 
proved computational effort. We have used this to study 
the 5 = 4 bound states of two types of Heisenberg 5=1 
Y-junction. 
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